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A hydro-dynamic description is used to study the normal modes of a vortex in a zero-temperature 
Bose-Einstein condensate. In the Thomas-Fermi (TF) limit, the circulating superfluid velocity far 
from the vortex core provides a small perturbation that splits the originally degenerate normal modes 
of a vortex-free condensate. The relative frequency shifts are small in all cases considered (they 
vanish for the lowest dipole mode with \m\ — 1), suggesting that the vortex is stable. The Bogoliubov 
equations serve to verify the existence of helical waves, similar to those of a vortex line in an 
unbounded weakly interacting Bose gas. In the large-condensate (small-core) limit, the condensate 
QQ ■ wave function reduces to that of a straight vortex in an unbounded condensate; the corresponding 

^\ ' Bogoliubov equations have no bound-state solutions that are uniform along the symmetry axis and 

0^ ■ decay exponentially far from the vortex core. 
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I. INTRODUCTION 



Recent experimental demonstrations of Bose-Einstein condensation in dilute low-temperature trapped alkali gases 
jlpH rapidly stimulated detailed theoretical Q-|tJ and experimental ^,|) studies of the lowest-lying normal modes. 
These latter experiments confirmed the essential applicability of the Bogoliubov approximation JlO| ] , which assumes 
that most of the particles remain in the self-consistent condensate for temperatures well below the Bose-Einstein 
, condensation temperature. Although Bogoliubov's original work described only a uniform system, the formalism was 
soon extended to a vortex line |Tl| , |T2| ] and other more general nonuniform configurations p3[ . 

In a harmonic trap with anisotropic frequencies u a (a — x, y, and z), an ideal-Bose-gas condensate would have 
characteristic dimensions d a — ^/fi/Muj a , where M is the atomic mass. As the number N of particles grows, however, 
the repulsive interactions become important, and the condensate expands from the noninteracting mean oscillator 
radius do = (d x d y d z ) 1 /' i to a larger mean radius Rq- For Rq 3> do, the kinetic energy of a particle in the condensate is 
negligible relative to the trap potential energy and the repulsive interparticle potential energy, leading to the Thomas- 
. Fermi (TF) approximation for the condensate density \L4 j. This limit occurs when the condensate number Nq is much 
larger than the ratio do /a, where a is the effective s-wave interatomic scattering length (typically Jig] do ~ a few /im, 
a rj a few nm, and the experimental value No > 10 6 amply satisfies the TF criterion that No 3> 10 3 ). In the TF 
J-^ ' regime, the expansion ratio is given by Ro/do w 15(A ? ' a/do) 1 / 5 3> 1 fhfll . 

^ . The normal modes of a nonuniform condensate can be described with either the pair of quantum-mechanical quasi- 
J^j ' particle amplitudes u and v that obey the coupled Bogoliubov equations |I^]l1^ , |16|| , or the equivalent hydrodynamic 
amplitudes for the fluctuations n' in the particle density and $' in the velocity potential p| ]l7|jL8[ . In the TF limit, 
the latter description simplifies considerably for the low-lying normal modes, allowing accurate predictions of lowest 
eigenfrequencies of an axisymmetric trapped condensate |||],[|. 

In contrast, the creation and detection of a vortex in a trapped Bose condensate remains relatively unexplored. 
The condensate density for a vortex line on the axis of an axisymmetric harmonic potential has been evaluated both 
numerically jl9],^0| and with the TF approximation Using the Bogoliubov equations in the TF limit, Sinha 

pl| studied the normal modes of a condensate with a vortex, showing that the presence of the vortex causes only a 
small shift of order (do/Ro) 2 C 1 for a subset of the eigenfrequencies. In addition, a detailed numerical study has 
been carried out for for the specific trap configuration of Ref. [Q for Na/do < 50 (for the distinct but related 
example of an annular condensate, see Ref. (^?J and references therein). Finally, Zambelli and Stringari |^8) have used 
sum rules to calculate the splitting of the lowest (quadrupole) normal modes caused by the presence of a vortex. 

The present work uses the equivalent hydrodynamic picture to study the low-lying normal modes of a quantized 
vortex line on the symmetry axis of an axisymmetric harmonic trap. In the TF limit, the vortex core radius £ ~ 
d^Ro is small compared to both the mean oscillator length do and the mean dimension Ro of the condensate. 
The resulting hydrodynamic equations then yield a perturbative treatment of the normal-mode eigenfrequencies and 
eigenfunctions for the density fluctuations, in which the vortex's circulating velocity splits the degenerate unperturbed 
eigenfrequencies with azimuthal angular momentum ±m by a relative amount of order £/i?o ~ (do/Ro) 2 *C 1, largely 
confirming Sinha's conclusion ]21fl based on the Bogoliubov equations (for \m\ = 1, however, we find additional terms 
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which ensure that the shift vanishes for the lowest dipole mode) . As emphasized by Zambelli and Stringari |2q] , this 
splitting provides a sensitive test of the existence of a vortex. 

Since the vortex-free trapped condensate is stable (both experimentally and theoretically), these results suggest 
that the corresponding trapped vortex is also stable at zero temperature. In the presence of weak dissipation, however, 
a slightly off-center vortex presumably will slowly spiral outward and eventually disappear, as noted by Hess p9| ] and 
by Packard and Sanders ||(]] for superfluid 4 He and by Rokhsar [^2| for trapped condensates. 

Section II reviews the basic hydrodynamic formalism, focusing on a condensate containing a q-fold quantized 
vortex in an axisymmetric trap. For the normal modes cx e lm ^, a general variational principle demonstrates that the 
frequencies are manifestly non-negative (and hence stable) for rn = and for m 2 > 4q 2 (Sec. V gives a more detailed 
treatment for the remaining cases, but only in the TF limit). For a large condensate (the TF limit), Sec. Ill makes use 
of the small parameter (,/Rq to perform a perturbation expansion of the hydrodynamic equations, yielding a splitting 
of the originally degenerate modes with ±m. In Sec. IV, the relative frequency shift of all low-lying normal modes 
is evaluated analytically for a spherical trap, whereas only a subset of the modes allows a closed-form expression 
for a strictly cylindrical trap and for a general axisymmetric but anisotropic trap. Section V contains additional 
considerations on the stability of a vortex line, based on the full Bogoliubov equations, including a demonstration 
that all normal modes with k = have non- negative eigenfrequencies in the TF limit (and thus are stable). 



II. BASIC FORMALISM 



Consider a trapped condensate characterized by the equilibrium condensate wave function ^ . It satisfies the Gross- 
Pitaevskii (GP) equation 0Jl| 

(T + H r + 5 |*| 2 -/i)f = 0, (1) 

where T = — 7i 2 V 2 /2M is the kinetic-energy operator, V tr is the external trap potential, g — Airah 2 /M is the effective 
interparticle interaction strength (expressed in terms of the s-wave scattering length a, here taken as positive), and 
/i is the chemical potential. It is convenient to write ^ = |\f | e lS , where no = l^l 2 is the static condensate particle 
density and vo = (7i/M)VS is the static condensate velocity. These quantities satisfy the static conservation law 
V-(novo) — 0, which is the imaginary part of the GP equation. Correspondingly, the real part 

pi = \^\~ x T\y\ + ±M V 2 + V tI + . 9 |vP| 2 , (2) 

generalizes Bernoulli's equation to include the quantum (kinetic-energy) contribution \^f\~ T\^\. 

Hydrodynamics relies on the local density and velocity as the relevant variables. In the present context, the linearized 
hydrodynamic equations determine the harmonic fluctuations in the particle density n'e~ tuJt and the velocity potential 
<I>'e~ Mt . The corresponding (complex) normal- mode amplitudes satisfy the coupled equations ||[[7j 

iwri = V- (v n') + V- (n V$') , (3a) 
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iu& = v -V$' + — n' %, — V- 

M 4M 2 n 



n V 



n 



(3b) 



both of which involve the (real) static condensate density no and (real) condensate velocity vo. These very general 
equations can be derived from either the linearized quantum-field equations or the conservation of particles and 
Bernoulli's equation for irrotational isentropic compressible flow Jl^Jl^]. The last term in Eq. (3b) corresponds to 



the kinetic energy (quantum) pressure and is omitted in the TF limit of a large condensate. Although the resulting 
formalism is frequently known simply as the "hydrodynamic approach," we prefer to use the term "hydrodynamic 
equations" for the more complete Eqs. @. 

The appropriate normalization of these amplitudes follows from that of the Bogoliubov amplitudes, which satisfy 
the condition flf] J dV (\u\ 2 - \vf) = 1 . A straightforward comparison of the linearized quantum-field equations in 
the presence of condensate flow jl7| shows that 

n = - (4a) 
*'=2MW^ + ^- ^ 

As a result, we have 
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= iM. („'*$' - $'V) 



(5) 



so that the normalization of the hydrodynamic amplitudes requires 

/ dVi(n'*& -$'V) = -^r. (6) 

Multiply Eq. ( |3a| ) by $'*, multiply Eq. (|^) by n'*, integrate the results over all space, and subtract the first from 
the second. A simple integration by parts yields 



<j dVi(n'*& -$'V) = J dV 



(n'*V -V& + n'v -V$'*) + n I V*'| 
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In this way, the normal-mode frequency is expressed as the ratio of two manifestly real integrals, so that u> itself must 
be real [assuming that the normalization integral Eq. (^) is nonzero] . In addition, variation of Eq. (Q) with respect to 
each of the two hydrodynamic variables nf* and <&'* reproduces the dynamical equations (||), so that Eq. (Q) serves 
as a variational principle for the normal-mode eigenfrequencies. Finally, Eq. (|^) shows that the left-hand side here is 
just Huj/M for a properly normalized set of solutions. 

We consider a quantized vortex with q quanta of circulation on the symmetry axis of an axisymmetric trap with 



V tt = \M{u\p 2 + uj 2 z z 2 ) = \Mu\{p 2 + AV), 



(8) 



where (p, <fi, z) are the usual cylindrical polar coordinates, and we define the anisotropy parameter A = uj z /lu±. The 



appropriate solution of the GP equation (P has the form * 
written v = V(j>, where 



V 



hq 
Hp' 



\^q{fi, z)\, so that the condensate velocity can be 



(9) 



In this case, the linearized hydrodynamic equations (|3a|) and (|3t]) have single-valued solutions of the form 

n'(r), <&'(r) oc exp(zm<^), (10) 
where m is an integer and the corresponding amplitudes depend only on p and z. Substitution into Eq. yields 
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where 



\^ q \ 2 is the (axisymmetric) condensate density for a g-fold vortex. 



Simple manipulation yields the equivalent form 
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(12) 



Here, the right-hand side is manifestly non-negative for m = and for m 2 > 4q 2 , so that all such hydrodynamic 
modes of a general g-fold quantized vortex necessarily have non-negative frequencies. As discussed in Sec. V, the 
quantum-mechanical "grand-canonical" hamiltonian in the Bogoliubov approximation reduces to a set of uncoupled 
harmonic oscillators for the different quasiparticle modes and is thus bounded from below so long as none of the 
eigenfrequencies is negative. Consequently, the vortex is stable with respect to modes with m — and m 2 > Aq 2 . 
For the remaining modes with < m 2 < Aq 2 , Eq. ( |l2| ) contains terms of both signs that can in principle produce an 
instability (in the usual case of a singly quantized vortex with \q\ = 1, this condition affects only modes with |m| = 1). 

Rokhsar has argued that a multiply quantized vortex is indeed unstable because of the formation of a bound 
state of radius ~ £ inside the vortex core (for a singly quantized vortex, he finds such a bound state only for relatively 
small condensates with Na/do < 5, where do is the mean oscillator length for the trap; a variational study of the 
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Bogoliubov equations confirms this qualitative conclusion). In the TF limit, the core radius £ becomes small and the 
details of the trap potential become irrelevant for phenomena on the scale of the vortex core, and the problem reduces 
to that of a long vortex in an unbounded condensate. For this case, the radial eigenfunctions for normal modes that 
are independent of z do not decay exponentially on the length scale £ p3|] , precluding such bound states as solutions 
of the relevant Bogoliubov equations (see Sec. V for a more detailed discussion). 

Although the variational Eq. ( |l2|) makes no use of the TF approximation, it does assume that the hydrodynamic 
variables provide the appropriate physical description. For a uniform condensate in the Bogoliubov approximation, 
the hydrodynamic picture holds only for long-wavelength modes with fc£ < 1 when the quasiparticle creation operator 
is an essentially equal admixture of a particle and a hole creation operator ^4j,^5|. For fc£ > 1, in contrast, the 
quasiparticle creation operator becomes a nearly pure particle creation operator. 

Correspondingly, a nonuniform trapped condensate has no hydrodynamic modes unless it is sufficiently large that 
£ < i?o, since the system otherwise approximates an ideal noninteracting gas. Even for such a large condensate, 
however, the only hydrodynamic modes are those with Tim < fj,, where p: is the chemical potential of the condensate 
j| . This latter condition becomes progressively more restrictive as A decreases and the condensate approaches the 
ideal-gas limit. 



III. PERTURBATION ANALYSIS OF THE NORMAL MODES 

We return to the hydrodynamic Eqs. (JsJ) . In the present case of modes cx e %m ^, they reduce to the coupled equations 



mV\ . m 2 , / 

n' - — + V- (n„ V$') = 0, 
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= o, 



(13a) 
(13b) 



where n„ is the condensate density for the g-fold quantized vortex. 



A. Stationary condensate 



Any practical solution of these equations requires a physically motivated approximation, and we here use the TF 
limit of a large condensate, when the kinetic energy in the GP equation (|j) is negligible compared to the remaining 
terms. In this limit, the equilibrium condensate density no = \^\ 2 is given by 

9 m 2 = n- 1 1 Mv 2 -V tr , (14) 

wherever the right-hand side is positive, and zero elsewhere; the chemical potential is determined by the normalization 
condition 

A w Nq — ( dV\^\ 2 (15) 



(the present zero-temperature approximation neglects the difference between the total number of particles N and the 
number Ao in the condensate). 

In the simplest case of an axisymmetric trap potential (||) and a stationary condensate with vo = Jl4| , Eq. jl^ ) 
gives the condensate density 

gn = uo - V tI , (16) 

where u is the chemical potential for a stationary condensate. The condensate density varies quadratically and can 
be rewritten in terms of the central density no (0) = no/g to give 



n (0) V R i R V V R ± R2 z 

where 0(x) is the unit positive step function. The condensate is ellipsoidal with radial and axial dimensions R± and 
R z given by 
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2 2pp 2^0 , 2 nS N 
Ra ~M^- h^^ (18) 



In the TF limit, each of these dimensions is much larger than the corresponding oscillator lengths d a = ^h/Muj a , so 
that fiQ ^> hiv a . Note that the aspect ratio is given by R z /R± = 1/A, with a cigar-shaped (disk-shaped) condensate 
corresponding to the limit A 1 (A ^> 1). 

The normalization condition (|lq) readily yields 



do 15 V hcjQ J 15 V do 



» 1, (19) 



where Rq = (i?^i? 2 ) 1 / 3 , d — (e^e^) 1 / 3 , and ujq = (^Wz) 1 / 3 are appropriate "geometric" means. It is conventional 
to define the coherence length £ in terms of the central density: 

^ 87mo(0)a ^ 



It satisfies the simple relation 

p 2 = 

2p M R 2 



,i K 2 di 

? = = 13> (21) 



implying the set of TF inequalities £ <C do <C Ro- 



B. Rotating condensate with a vortex 



In the presence of a quantized vortex with q quanta of circulation on the trap's symmetry axis, Eq. (jlj) yields the 
corresponding condensate density 



gn q = H q — V tI ~ \MV 2 , 

where p q is the chemical potential for a q-fold vortex and V = fiq/M p. This expression can be rewritten as 

n q (p,z) - 2 ~ 2 » 2 ~ 2 



n (0) 



Pg P 
Po R± 



z- q 2 ^ 



c( _ fL _ il _ i!i 2 



(22) 



(23) 



22|. The normalization 
P q - Po)/po for fixed 



using Eq. (|2l|). The centrifugal barrier arising from the vortex alters the shape from ellipsoidal to toroidal; in the 
TF approximation, it creates a flared hollow core with radius ~ \q\£ in the plane z = [^o] 
condition ( |l5| ) can be expanded to show that the fractional change in the chemical potentia 
V tI and N is of order q 2 d^/R^ In (R 2 /\q\ df) < 1 @- 

Equation (|l3b|) ca n be simplified considerably in the TF limit when the kinetic energy is treated as small ||. 
Specifically, Eq. ( |13q ) contains the density fluctuation n' in two physically distinct ways. The last two terms involving 
h 2 /M 2 n q reflect the quantum kinetic-energy contribution to Bernoulli's equation; in contrast, the preceding one 
involving g/M reflects the repulsive interparticle interaction energy. In the present case, the kinetic-energy term for 
p > £ is smaller than the interparticle potential energy term by a factor of order ft 2 /gMrio(0)R 2 « £ 2 /R\ 



As will be seen below, the effect of interest here is of order £ / R 
can be omitted entirely for p > £, giving 



d^/R 2 , so that the kinetic-energy term in Eq. (p_3q) 



mhq 



0. 



Mp 2 ) ' M ' 

A combination of Eqs. (|l3a|) and (E4) yields the single eigenvalue equation 



gn q 

M 



M 



mhq V 
Mp 2 J 



= 



(24) 



(25) 



This second-order partial differential equation determines the frequencies and normal modes of an axisymmetric 
trapped condensate containing a vortex line through terms of order d\fR\. A combination of Eqs. ([l^) and ( |l9| ) 
shows that this small parameter has the following dependence 
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15iVaA J 



(26) 



on the condensate number N and the asymmetry parameter A 

It is now convenient to introduce dimcnsionless units, with u>' = uj/lu±, p' = p/Rj_, and z 
in particular, the dimensionless ratio rnhq/tuMp 2 = (mqdl/ R^p'^ 2 is seen to be of order d 2 
In addition, the actual condensate density n q differs from that for a vortex-free condensate 



= z/R z = zX/Rj 



JR 
n 



j_ and thus small, 
by corrections of 



order £ /Rq = d^/R^, which is negligible compared to the effect of the superfluid flow. Omitting the primes on the 
dimensionless variables, we therefore find 
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2n (0) 



ld_ 

pdp 



pn 



9$' 
dp 



m 2 n <&' ^ 2 & 
p 2 dz 



n - 



~dz~ 



2ujmq d 2 
~p 2 ~R z 



$' = 0, 



(27) 



where terms of order d^_/R\_ have been omitted, and the TF condensate density no is taken from Eq. (|l7j). 

Equation for the velocity potential <&' contains an explicit term of order cij_/i?j_ arising from the circulating 
condensate velocity. In fact, however, we seek solutions in a class of functions for which the particle-density fluctuation 
n' is finite on the vortex axis (at p = 0). To incorporate this boundary condition, it is natural rewrite Eq. (|27 

terms of the function n'. Equation ( |2~i| ) shows that n' cx fl — Z7pl?~) < ^ > '' or ec l u i va l en tly, 



$' cx 1 - 



mqd\ 



LUp 



mqd\ 
2 R 2 



up 



Then, omitting terms of order d\/R\, we can replace Eq. d27j) for $' by the following equation for n': 
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drt_ 
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p 2 dz \° dz 



+ LU 2 n' + 
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«o(0)/0 \p 



d_ 

dp 



= o 



(28) 



where no/no(0) = 1 — p 2 — z 2 in dimensionless units. 

The normal modes of a large vortex-free (q = 0) axisymmetric condensate satisfy the following equation ]3l|j3S 
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2n (0) 
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dn' n 



(29) 



where up is the frequency of a particular unperturbed normal mode. These low-lying TF normal modes have been 
studied analytically for a spherical Q and a cylindrical trap |33j ; for a general axisymmetric trap, some of them can 
be found analytically p,pl| , p2| , but, in general, numerical techniques are necessary Q. 

Comparison of Eqs. (|28|) and ( p9| ) immediately yields a simple perturbation expression for the frequency shift of 
any particular normal mode of an axisymmetric condensate induced by a g-fold quantized vortex 



or 



2mqd\ I (w ) 2 

urn 
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n Q {0)p 3 



d_ 

d~ P 



(30) 



where the angular bracket denotes a matrix element evaluated with the appropriate eigenfunction of the vortex-free 
condensate. Note that this result describes only modes with m/0 because the convective operator v • V = (V/p)d/d(j) 
vanishes for an axisymmetric state with m = 0. To the accuracy considered here, the region far from the vortex core 
makes the main contribution to the frequency shift, and the result is independent of the detailed core structure. 

One should mention that if we take Eq. @ for instead of Eq. (fH for n', as a basis for perturbation analysis, 
we obtain equivalent results for the frequency shift only for |m| > 2. The point is that the unperturbed solutions for 
$ and n' both behave like p' m ' near the vortex line and, as follows from Eq. (HJ), n' will finite at p = for such 
solutions only if |m| > 2. But if we treat Eq. (^8|) for the particle-density fluctuation n' perturbatively, the boundary 
condition for n' is satisfied automatically for any value of m. As a specific example that illustrates this distinction, 
Eqs. (58) and (|59|) show that the exact lowest dipole modes of a vortex (those with |m| = 1) have n' cx p and <&' oc p^ 1 
for p — > 0. 
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IV. FREQUENCY SHIFTS FOR SPECIFIC CASES 



It is convenient to write the fractional shift in the squared frequency as 

,2 _ (,.,0\ 2 J2 



where sgnm = m/\m\ and 



or — lor i . d, , , 

-^^Asgnr,^, (31) 



A (J°)3 ( ( p 2 + n (0)p 3 (dp pj ' ' 



Since the matrix element is independent of the sign of to, the vortex, generally speaking, splits the two previously 
degenerate modes with ±|to| into two levels, according to the sign of to. In particular, the (small) splitting for \m\ > 
can be characterized by the expression 



\m\ ~ W-M _ A d ± A ( d± \ 



- = ?A [ T71 ^ rT ) , (33) 



where the last form follows from Eq. (pq). The remaining calculation depends on the shape of the harmonic trap, and 
we consider several specific cases. 

A. Spherical trap 

For a spherical trap, each unperturbed eigenfunction is a product of a spherical harmonic and a radial polynomial 

<oc^ m (0,0)r / P r ( l i+i '° ) (l-2r 2 ), 

where Pn +2 '°\x) are Jacobi polynomials, r 2 = p 2 + z 2 , p = rsin#, and n is the radial quantum number. The 
corresponding dimensionless unperturbed eigenfrequency is given by || (^ni) — I + n(2n + 21 + 3). 

To calculate the matrix element, one performs an angular and radial integration. The angular averaging of the 
corresponding operators gives rise to the following result 



sin 



1 \ 21 + 1 



2\m\ 



2 , 



(34) 



<2 



\ m]>1 , 



COS0 d 1 \ _ I 4|m| 

^ede~^~e/ n ~y _(2i+il (l + ffi+i)), | m | = 1 . 
The remaining radial integral can be evaluated using the following integrals for Jacobi polynomials |}5j 

1 (1 - xT[P^\x)Y = - (36) 

_i a + 2n + 1 

(l-z)«-^ a ' 0) (s)] 2 = ? (37) 



L (1 _ ,)^[P(-0) (X)]2 = 2^(a 2 +2na + a + 2n 2 + 2n) 
i ' a (ar — 1) 

As a result, the general frequency shift for the spherical mode with quantum numbers (nlm) is given by the expression 
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f . 4 " +2 '+ 3 ■, |m|>l 

A nZm — \ 4»+2Z+3 _ (4n+2i+3)Z(7+l) [(1+2tQ (i+3/2)+2n 2 ] 

I % /i+n(2n+2i+3) _ 4(7-1/2) (1+3/2) [l+n(2n+2l+3)] 3 ' 2 

Equations (|3^) and ( p9| ) show that the vortex splits some (but not all) of the previously degenerate 21 + 1 modes with 
given (nl) and different m. The m = mode remains unshifted in this approximation, the |m| = 1 mode is split by 
an amount A^ x , and the remaining modes with 2 < \m\ < I are split by a different amount A^ m that is independent 
of the particular value of m. For |m| = I = 1, n = (the lowest dipole mode), one can see that the frequency shift 
and splitting vanish, as expected, because the vortex-free condensate and that with the vortex both oscillate at the 
bare trap frequency loj_. 



(39) 



B. Cylindrical trap 

The cylindrical case is obtained from the general trap potential by setting ui z = 0, but it is also necessary to use 
the radius Rj_ as the characteristic length for both the radial and axial coordinates. The unperturbed condensate 
density is now given simply by no(p) = no(0)(l — p 2 ), and the perturbed eigenfunctions can be written in the form 
n' {p)e % ^ m ^ +kz \ where n' obeys the eigenvalue equation 



1/1 2^ f l d d m ,2n 



dp 



n'(p) 



2mqd 2 
ujp 2 R 2 



'2 r i ^2 



2 , ,-p 2 fl d 



p dp 



n'(p) = 0. (40) 



The term with l^R 2 ^ precludes an exact solution [ p3[ , but the frequency shift can again be found perturbatively for 
wavelengths long compared to R± , when we find 

2 . ,2_ 2mqd\ /(cj°) 2 -1 , 1-p 2 ( d 1\\ 12 



•"^v^+vni-^+i*^ 1 -^- (41) 

Here, the angular brackets denote a matrix element evaluated with the unperturbed eigenfunctions, which have 
the explicit form n' cx p\ m \p^ m \'°\l — 2p 2 ), and the corresponding eigenfrequency is given by [j31p5| l (w° m ) = 
\m\ + 2n(n+\m\ + l). 

The fractional shift in the squared eigcnfrcqucncics is readily evaluated with the same radial integrals used for the 
spherical case 



to 



2 - (W° m ) 2 k 2 R\ d\ ^\rn\+2n(n+\m\ + l) 



(4 " +2 l m ' +2) |m|>l 



2(\m\ + 2n)(\m\ + 2n + 2) q& R 2 ± ) 4n(» + i)(n+ 2 ) , | = 1 

[l+2n(n+2)] 3/2 ' 1 1 

where the second term is the contribution of the vortex and the first is that of the traveling wave. The vortex again 
has only a small effect, which implies that the configuration is stable, assuming that perturbation theory is valid. For 
\m\ = 1, n = 0, k = (dipole mode), the frequency shift is equal to zero. 



C. General axisymmetric trap 

The normal modes of a general axisymmetric trap in the TF limit have been classified completely and reduced to 
the diagonalization of a sequence of increasingly large matrices [^l|,^2| . We shall use the cylindrical coordinates of 
Rcf. j3^|, which turn out to be more convenient for the present purpose. The unnormalized solutions of Eq. (^9|) have 
the form n'^~(p, z) — p' m 'B^ m and n'^~ (p, z) = p\ m \zB~ m and are, respectively, even and odd under the transformation 
z — > — z (our notation differs somewhat from that of Ref. p^|). Here, B + m and B~ m are even polynomials with 
typical terms p 2ni z 2n2 , where ni, ri2, and n are non-negative integers with n\ + n,2 < n. For any given n and m, 
there are n + 1 even-parity normal modes and, separately, n + 1 odd-parity normal modes, labeled by j = 0, 1, • • • , n, 
with frequencies determined by diagonalizing (n + l)-dimensional matrices given in Ref. []32fl . For n — 0, the 

even (+) and odd (— ) modes have dimensionless frequencies uJ®Q m — y/\m\ and Lu®Q m — y/\m\ + A 2 , respectively, 
with Bq 7H = B^ m = 1 and j = 0. For general values of the asymmetry parameter A = oj z /oj±, it is straightforward 
to determine the corresponding solutions for n — 1 p2| , but the higher-order modes require numerical analysis (for 
A —>■ 0, Rcf. shows that all the unperturbed frequencies agree with those for a cylindrical trap). 
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Equations ( J3l| ) and ( |32| ) determine the vortex-induced frequency shift A- nm of any particular normal mode. The 
matrix clement should be evaluated for j = 0, 1. 2, ■ ■ ■ , n (the n + 1 distinct eigenfunctions for a given nm±). For the 
simplest case of n — 0, we find 

2|m|+3 I I . -t 
! — m > 1 



A + 0m =<| ^ H ' (43) 

), \m\ = 1 



2|m|+5 I I . 
m > 1 



A = { ( 44 ) 

(i+A 2 ) 3 / 2 ' |m| — 1 . 



These expressions are equivalent to those found in Ref. |21| apart from the case \m\ = 1, where our treatment of the 
singular terms gives additional terms, ensuring that the lowest dipole mode is unshifted because A^ = 0. 

As a check, we note that Aj 0m is just from Eq. ( |39| ) for a spherical condensate with n = and / = |m|, 

and that Aqq to for A = 1 is just A^ m for a spherical condensate with n = and I = \m\ + 1. In addition, Eq. ( |3^ ) 
with our values for Aq 02 = 7/V% and Aq 01 = 7A 2 /(1 + A 2 ) 3 / 2 yields frequency splittings that agree with those found 
independently from sum rules by Zambelli and Stringari ]28| for the lowest quadrupole modes with \m\ =2 and 
|m| = 1 in the TF limit. 

For the next case of n = 1, Figs. |l| and [2] show the corresponding dimensionless eigenfrequencies and fractional 
frequency shifts and splittings A^ lm for the even and odd modes j = 0,1 and \m\ — 1, 2 as a function of the asymmetry 

parameter A at fixed d\/R\. It is interesting that the curves for A^ xl and are non-monotonic functions of A 

and reach a maximum near A = 1, which is a spherical trap (this behavior is also seen in Aq 01 above). In contrast 
to the situation for the unperturbed frequencies, the frequency shifts A* for a long cigar-shaped axisymmetric trap 
(the limit A — > 0) differ from those for a cylindrical trap because of the different angular integrations. To check these 
results, the particular case of A = 1 (a spherical trap) was studied analytically for n = 1, confirming that these normal 
modes are indeed consistent with those found in Ref. 01. 



V. BOGOLIUBOV EQUATIONS AND VORTEX STABILITY 

In the Bogoliubov approximation, the "grand-canonical hamiltonian" operator K = H — [iN reduces to a pure 
condensate part Ko that depends only on the condensate wave function \t, and a noncondensate part K' Jl^Jl7[| . The 
latter involves the solutions of the Bogoliubov equations, 



Ca, ( /v|'-Y„ /•;,//,. (45) 



Cvj-gi^fu^-E^, (46) 



where C = T + Vt r + 2g|5'| 2 — ji. Here, Uj and Vj are a complete set of coupled amplitudes that obey the normalization 
condition J dV {\uj\ 2 — \vj\ 2 ) — 1, and Ej is the associated eigenvalue. The ground-state occupation of the jth mode 
is simply iV' = J dV\vj\ 2 , and the noncondensate part (a second-quantized operator) reduces to 

i'=-^£ 3 iV;+^Vjft, (47) 

J 3 

where the primed sum means to omit the condensate mode, and /3j and (3j are quasiparticle creation and annihilation 
operators that obey boson commutation relations [f3j,f3 k ] = Sjk- The ground state of K' is the vacuum |0), defined 
by /3j\0) — for all j ^ 0, so that the first term in Eq. ( |47| ) is the ground-state expectation value of K'; the remaining 
term involves the quasiparticle number operator (3j(3j, and the spectrum of K' is bounded from below as long as 
Ej > for all j. If, however, any of the properly normalized solutions has a negative eigenvalue, the system can 
lower its ground-state thermodynamic potential Q = (K) arbitrarily by creating more and more quasiparticles for 
that particular normal mode. This behavior indicates that the original ground state 4" is unstable and must be 
reconstructed to form a new stable ground state pi] . 
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A. Axial vortex waves 



For a condensate containing a q-iold quantized vortex in an axisymmetric trap (studied in Sec. II), the condensate 
wave function has the form ^(r) = e lq<t '\'$> q (p, z)\. In the presence of a general condensate velocity v = (7i/A/)VS, it 
is convenient to make explicit the phase of the condensate wave function through factors e ±lS jlTj , and the Bogoliubov 
amplitudes here can be written as 



u ( r ) \ _ imip+ikz 

v(r) 1 ~ ' 



e-***v m (p,z) 



(48) 



This form represents a helical wave with wavenumber k and angular momentum to relative to the condensate; we 
assume that the factor e zkz varies rapidly compared to the amplitudes u and v, with the condensate containing many 
wavelengths (kR z ^> 1). The corresponding amplitudes obey the coupled equations |l2j 



2M 
2M 



V 2 + 



2M 
2M 



(to + q) 2 
P 2 

(to — q) 2 



+ k 2 



2ik 7T 
oz 

2ik 7T 
oz 



V t r + 2g\V q 



Vtr + 2g\$ q \ 2 - p 



\ 2 v, 



g\* q \ 2 ur, 



Eu r , 



= —Ev„ 



(49a) 
(49b) 



where V 2 = d 2 /dp 2 + p 1 d/dp+d 2 /dz 2 is expressed in cylindrical polar coordinates. The different centrifugal barriers 
imply that the two amplitudes behave differently near the axis of symmetry, with u m oc and v m oc pl m ~ ? l as 

As noted in Sec. II, these Bogoliubov amplitudes are linearly related to the corresponding hydrodynamic variables 
n' and In an axisymmetric trap with a g-fold vortex on the symmetry axis, their angular dependence is simply 
e"™*, and the corresponding amplitudes are given by Eqs. (|j) 0,^8 37 

n 'jm = 1*9 1 («jm ~ , (50a) 



^ m = 2Mi\^! q \ + ' (50b) 

where j denotes the remaining set of quantum numbers, and |^q| cx p' 9 ' for p — > 0. In the simplest case of a vortex- 
free condensate (q = 0), these amplitudes reproduce the behavior cx p'™ 1 ' found in Refs. |3^J32], but the situation 
is more complicated for q =/= 0. The boundary condition that u and v remain bounded near the origin (discussed 
below in Sec. V.B and in Ref. p3| ]) implies a corresponding behavior for n' and Although $' has an apparent 
singularity for small p from the factor l^gl -1 , physical quantities like the hydrodynamic current involve an additional 
factor n q = \^ q \ 2 p7| . As shown in Sec. Ill, the detailed core structure has negligible effect on the normal-mode 
eigenfrequencies in the TF limit (a large condensate), but it would be significant for a smaller condensate with 
Na/do ;$ 10, when the TF relation in Eq. (19]) implies that do/Ro ^ 0.4 is no longer small | pr| . 

The general Bogoliubov equations are difficult to solve. In the present case of a harmonic trap potential Vt r , 
however, they have three exact solutions. For these special "dipole" modes, any condensate described by a wave 
function \E' that satisfies the GP equation (|l|) undergoes a rigid oscillation along the three principal axes of the trap 
with the bare oscillator frequencies |[7[ . Dehne the raising and lowering operators 



1 f x a d 

a a = T" + da "5 

a/2 V "a OX a 



(51b) 



where d a = ^Jh/muj a , and a = x, y, and z. As shown in Ref. |37|, the state 



v a (r) J \a a V*(r) 
satisfies the coupled Bogoliubov equations with an eigenvalue E a — huj. 



U*(r)^iy£l) = (f a Z { rl) (52) 
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For an axisymmetric trap potential Vt r = iM w 2 p 2 + \Muj\z 2 , the two states U x and Wj, are degenerate, each with 
energy hu>±. In terms of the new operators |p8| 



a± = — (a x =F ia y ) , 
a ] ± = — (4 ± «4) ' 



(53) 
(54) 



the (conserved) angular momentum about the symmetry axis can be rewritten as L z = h(a? + a+ — a_a~). As a result, 

(a\_ ) adds a unit of positive (negative) angular momentum. Correspondingly, we have two distinct solutions of the 
Bogoliubov equations 



W±(r)= ( , ) =U x (r)±iU y (r) 



<4*(r) 
a T **(r) 



(55) 



each with energy ?kj_L. In cylindrical polar coordinates, the explicit form of these amplitudes is readily determined; 
for a g-fold quantized vortex, we find 



u + (r) 
77+ (r) 



2 e 



2 e 



P , d qd A 

i — d± ir + — 

d± op p 

>, d qd ± 

- + d± — H 

1 dp p 



l*«(p.*)l, 



(56a) 
(56b) 



U-(r) 
U-(r) 



l e i(-3-l)* [ JL 



9 gdj 
"dp p 



I* 9 (P. 



(57a) 
(57b) 



and it is clear that these solutions are properly orthogonal, with J dV (u*_u+ — v*_ W+) = 0. Comparison with Eq. ( J4S 
shows that W± has m = ±1 and fc = 0. Pitaevskii |H| constructed analogous zero-frequency solutions in the case of 
a vortex line in an unbounded fluid (his solutions emerge in the limit wj_ — > 0, so that d± — * oo). 

It is worth noting that the fluctuations in the particle density n' and the velocity potential $' for these lowest 
dipole-mode solutions (p6h-(E^) have a simple explicit expressions: 



n = e 

2 



dn q 
dp' 



(58) 



2Mi 



P_ qd±_ 

d±_ p 



(59) 



The contrasting behavior of n' and $' near the vortex axis helps motivate the discussion in Sec. III.B. 

The form of Eq. (49) suggests an approximation for small k and m = ±1, using U± as the unperturbed eigenfunc- 
tions. Standard first-order perturbation theory gives the corrected eigenvalue 



E± (k) « ?ia;_L + 



fr 2 fc 2 J dV (\u ± \ 2 + \v ± \< 
2M JdV(\u±\ 2 - |u±| 5 



(60) 



where the additional term involving d/dz yields a surface contribution that vanishes. Use of Eqs. ( p6|) and ([57]) yields 
the explicit expressions 



/ dV(\u ± \ 2 + \v ± \ 2 ) = 2 y 



P , d± 
— ±q — 

dj_ p 



l* 



9p 



(61a) 



j dV(\u ± \ 2 ~\v ± \ 2 )=2 J dV (p±<Zy 



9 P 



(61b) 



ii 



Equation (61b) is easily evaluated with an integration by parts and the normal izati on condition Eq. ([15]). In the TF 
limit of a large condensate with d± <C R±, the dominant contribution to Eq. ( 61a ) comes far from the vortex core, 
so that we can replace \H! q \ from Eq. ( p3| ) by \^ \ from Eq. (17). Neglecting corrections of order d^/R 2 ^, we find the 
compact result 



E±(k) « hujj_ 



h 2 k 2 

2M 



^±,)=^[1 



\ {\k 2 R 2 L ±qk 2 d 2 L ) 



(62) 



As in Eq. (|42|), the correction is indeed small if \q\ is not too large and kR± < 1 (since we also require kR z ^> 1, this 
latter condition implies a cigar-shaped condensate with small asymmetry parameter A <C 1). Consequently, a q-fold 
quantized vortex has long-wavelength propagating axial helical modes with wavenumber k and angular distortion 
m = ±1. In this perturbative limit, the frequencies of these two modes exceed the bare dipole frequency uj±, as 
expected for an internal distortion. The separate parabolic branches are split by an amount 



E+(k) - E_{k) » q 



h 2 k 2 

M ' 



(63) 



In his original work on a singly quantized vortex in an unbounded dilute Bose gas, Pitaevskii jl2j developed an 
approximate description of vortex waves for shorter wavelengths, assuming only that fc^<Cl. To logarithmic accuracy, 
the dispersion relation has the classical form uik = (hk 2 /2M^ In (l/fc£) [^9|. This method is easily generalized to the 
present case of a g-fold quantized vortex on the symmetry axis of an axisymmetric trap in the TF limit. We use the 
coherence length £ as the length scale and rewrite the exact condensate density as ip 2 = \^ q \ 2 /no(0) . The coupled 
Bogoliubov equations (E9|) for m = — 1 become (here, the coordinates are dimensionless) 



oz a 1 H 



-V 1 



P 

(q + l? 



U-l — ij} q v^i = CU- 



V-l - tpnU-l 



-ev_i, 



(64a) 
(64b) 



where we have set fi q f« = 5^o(0) and introduced the abbreviations e = (2E /Tilo^_)(^ 2 / d\) = E/no ^ 1 and 
k = fc£ < 1. 

The amplitudes vary slowly along the axis with a characteristic length i? z , so that the various derivatives with 
respect to z can be omitted, along with the trap potential (which is here of order £ 4 /c^ <C 1), in which case these 
equations (64) take the approximate form 



1 d d (q-lf 
p dp dp 



1 d d 
pdp^ dp 



(q + l? 
P 1 



2< - 1 



2^-1 



- ip„U-x = -ev- 



(65a) 
(65b) 



In the same approximation, the condensate wave function satisfies the following GP equation 

1 d d q 2 , 2 \ , 

-pdp-% + 7 + ^' 1 ) i;q = > 



(66) 



whose asymptotic solution for p 3> 1 is tp q f« l — q 2 /2p 2 . This set of equations is now precisely the same as that obtained 
by Pitaevskii |]l2| , apart from the appearance of the quantum number q instead of q = 1. Exactly the same arguments 
then show that the dispersion relation is given by E^ w q(h 2 k 2 /2M) ln(l/fc£), where we assume t;^Cl< kd± . This 
condition implies that Ek is significantly greater than the lowest dipole energy E± — huj± = Ti 2 /Md\. 



B. Bound core states 



Rokhsar p2| has argued that the mere existence of a bound state localized in the vortex core implies an instability, 
with the original axisymmetric configuration undergoing a displacive transition of the vortex core. In the present 
context, such an instability would arise only from the presence of one or more negative-energy eigenvalues of the 
Bogoliubov equations, causing the original ground state of Eq. (47) to collapse. For such a situation, our preceding 
perturbative approach would fail completely, requiring a more powerful approach. 
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In their original solution for a long vortex line in unbounded superfluid 4 He, Ginzburg and Pitaevskii ]40| remarked 
that multiply quantized vortices are energetically unfavorable. It is also notable that no clear evidence for doubly 
quantized vortices in superfluid 4 He II has ever been reported [ ElJ . Thus, we consider only the case q — 1 and seek 
solutions with relative angular symmetry oc e"™^ for the full Bogoliubov equations with a singly quantized vortex in 
an axisymmetric trap. In the TF limit, we neglect small corrections of order £, 2 /Rq ~ £ 4 /^o> and these equations 
then reduce to those for a vortex line in an unbounded condensate [here, they are direct generalizations of Eqs. (|6^), 
obtained for general m by setting q = 1]. For q = 1, the result of Eq. ( |l~2] ) shows that the only possible unstable 
solution corresponds to the choice |m| = 1. Then, since the amplitude u-\ remains finite at the origin, we anticipate 
that the choice m = — 1 will yield a positive normalization for |u_i| 2 — | 1 1 2 - 

These coupled Bogoliubov radial equations ( |65|) are equivalent to a single fourth-order ordinary differential equation. 
Near the origin, there are four linearly independent solutions, and it is straightforward to show |23j ] that only two are 

finite as p — > 0. In a matrix notation, they can be written (u^l,v^l) and (u_l,V^l). Each is a power series in p 2 , 
with leading behavior (omitting constant factors) 

£8)-(Z) - 

where we have suppressed the subscript — 1. The general solution is a linear combination, with u « ci + 0(p 2 ) 
and i; « C2/3 2 + C(/9 4 ), where Ci and C2 are constants. Note that the linear combinations u ± v appearing in the 
hydrodynamic amplitudes [Eqs. (|50|)] both have the same behavior at the origin. Although this situation apparently 
produces a singularity in the variational expression in Eq. (|i~2|), it is not hard to verify that the divergence in fact 
cancels out. 



Similarly, it is not difficult to verify |23| that the asymptotic solutions of Eqs. ( ^5h have the form cx p^ 1 / 2 exp(i£p), 
where C 2 is the solution of a quadratic equation. One solution is negative with £ 2 = — Q m = — Ve 2 + 1 — k 2 — 1; 
the other is £f = \/e 2 + 1 — k 2 — 1, which can take either sign. To eliminate the growing exponential solution 
cx p -1 / 2 exp(|Ci m |p), it is necessary to adjust the ratio C1/C2, leaving a solution for u and w that is bounded both at 
the origin and at infinity. It has the asymptotic form 

u, v ~ p- 1/2 [Aexp«i / 9) +Bexp(-iCip) + Cexp(-|Cim|p)] , (68) 
where A, B, and C depend on e and k. 



For k 2 > 0, it is easy to see that £1 is real if |e| > €b{k) = Vk 4 + 2k 2 , which is the Bogoliubov energy for the 
wave number k along the vortex axis. For this range of parameters, it follows that the asymptotic solutions oscillate 
and there is no bound state. In contrast, if |e| < £b(k), then £1 becomes imaginary, and the coefficient of the 
growing exponential must vanish. If we choose Im £1 > 0, then B(n, e) must vanish, which in principle determines 
the dispersion relation of the bound state as a function of the axial wave number. In effect, this procedure yields 
Pitaevskii's vortex-wave solution |12| . 

The situation is different if k = 0, for then e^(0) vanishes, and £ 2 = \/e 2 + l — 1 is non-negative for all real e. Thus 
the general solution for k = necessarily oscillates far from the vortex core for any real positive value of e 2 , implying 
that there is no bound state in this case (the solution with e = is simply the dipole oscillation, here with zero 
frequency, as found by Pitaevskii p"2|). A similar argument holds in the extreme TF limit for general q and m, since 
there are always two linearly independent solutions that remain finite at the origin [p3| , and the asymptotic behavior 
for large p is independent of the values of q and m (the centrifugal barriers contribute only in higher orders) . 

This behavior is easily understood from a classical perspective. A long classical hollow-core vortex line in an 
incompressible fluid has irrotational flow that is describable with a velocity potential. In addition to boundary 
conditions at the inner surface of the core, the velocity potential satisfies Laplace's equation [ff2| . Solutions of 
Laplace's equation cannot exhibit spatial oscillations in all directions, so that a vortex wave propagating along the 
axis with wavenumber k necessarily decays exponentially ~ e~ kp in the radial direction. In the limit k — > 0, this 
exponential decay length k^ 1 diverges, and the amplitude can only vanish algebraically as p — > 00, which is the 
behavior found here for a vortex line in an unbounded condensate. 

Based on Rokhsar's discussion p3], the numerical study of Dodd et al. pfl found one negative-frequency normal- 
mode solution for a singly quantized vortex (the actual negative frequency should be the negative of that shown in 
their Fig. 2). In the ideal-gas limit, the anomalous frequency is — ui±, and it increases toward zero from below in the 
large-condensate limit Noa/do 1. The present treatment implies that this particular eigenfrequency approaches 
zero from below in the limit of a large condensate (£/i?o -C 1). The precise dependence on the small parameter, and 
the relation to the critical angular velocity for vortex creation remains for future analysis. 
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VI. DISCUSSION 



We have used a variety of methods to study the stability of an axisymmetric vortex in a large axisymmetric trapped 
Bose condensate. In this limit, the repulsive interparticle condensate energy dominates the condensate kinetic energy. 
The resulting Thomas-Fermi (TF) approximation holds when the equatorial condensate radius R± is much larger than 
the corresponding oscillator length d± = y/h/Mtj±. The hydrodynamic description of the condensate's normal modes 
simplifies greatly in the TF limit, and the presence of vortex leads to small corrections of relative order d\/ R\ <C 1. 
A perturbation analysis shows that a vortex splits the degenerate normal modes of the vortex-free condensate (those 
for ±m) by an amount of this same order, which should provide a clear signal for the existence of a vortex in the 
trapped condensate. Since the vortex-free condensate is considered stable from both a theoretical and an experimental 
perspective this conclusion indicates that the vortex should also be stable, at least in perturbation theory. To 

preclude the possibility of nonperturbative effects associated with the possible presence of bound states with a radial 
range ~ £ localized in the vortex core [z| , we demonstrate that a long vortex line in an unbounded condensate (the 
extreme TF limit of a trapped condensate) has no such normal modes p3[ ] other than the vortex waves found by 
Pitaevskii |||. 
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FIG. 1. The A-dependence (asymmetry parameter) of the unperturbed dimensionless frequency u°jt m and the fractional 
frequency shift A^ lm for even modes with n = 1 and j = 0, 1 (note that A < 1 represents a cigar-shaped condensate). Solid 
lines denote m — 1 and dashed lines denote m = 2. 



FIG. 2. The A-dependence (asymmetry parameter) of the unperturbed dimensionless frequency w° lm and the fractional 
frequency shift A~ lm for odd modes with n = 1 and j = 0, 1 (note that A <C 1 represents a cigar-shaped condensate). Solid 
lines denote m — 1 and dashed lines denote m = 2. 
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